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On Implicit Runge-Kutta Methods for 
Parallel Computations 

Stephen L. Keeling* 


Abstract. Implicit Runge-Kutta methods which are well-suited for parallel computations are 
characterized. It is claimed that such methods are first of all, those for which the associated 
rational approximation to the exponential has distinct poles, and these are called multiply implicit 
(MIRK) methods. Also, because of the so-called order reduction phenomenon, there is reason to 
require that these poles be real. Then, it is proved that a necessary condition for a qr-stage, real 
MIRK to be A-stable with maximal order q + 1 is that q = 1,2, 3, or 5. Nevertheless, it is shown 
that for every positive integer q, there exists a 9 -stage, real MIRK which is Ao-stable with order 
q + 1, and for every even q, there is a 9 -stage, real MIRK which is /-stable with order q. Finally, 
some useful examples of algebraically stable MIRK’s are given. 


‘Supported by the National Aeronautics and Space Administration under NASA Contract No. NAS1-18107 while 
in residence at the Institute for Computer Applications in Science and Engineering (ICASE), NASA Langley Research 
Center, Hampton, VA 23665-5225. 
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1 Introduction. 


This paper is concerned with the characterization and construction of implicit Runge-Kutta 
methods (IRKM’s) which are especially well-suited for the approximate solution of evolution equa- 
tions on a parallel machine with a modest number of processors. For a precise discussion of the 
issues, IRKM’s and their properties are now introduced. Given an integer q > 1, a q-atage IRKM 
is determined by a set of constants: 


an 

aj. 

n 

a 9 i 

‘ " a 99 

T i 

bi 

••• b q 



and it is convenient to make the following definitions: 


A = K}l=i> b T = (b u b 2 ,...,b q ), 



v - ' T i Uj=i> 


dig ft}, 

R = diag {i -1 }, 

1<»<9 


M=BA + A t B - bb T , 

e T = (l, 


For the IRKM formulation used in this work, choose arbitrarily, to G U, yo € R", F : R n+1 — > R n 
sufficiently smooth, and A: > 0 sufficiently small, so that for to < t < to + k, smooth functions 
y,y : R — + R n are well-defined by: 


I D t y(t) = F(t,y(t)) 

\ y{to) = yo, 

9 

yo + (* - *o)XX F ( f ° + t A l ~ to),y J (t ))> i <j<q 

;=i 

9 

yo + (i - fo£&iF(to + T A l ~ *o),y’(*))- 

»=i 

The method is described as explicit if aij = 0, i < j and implicit if for any t, a,-,- ^ 0. 

As usual, there are three criteria by which a method is judged in this work: order of consistency, 
stability, and implementability. First, an IRKM is said to have order v if for every y and y defined as 
above, D l t y(to) = Djy(to), 0 < / < v. Nevertheless, when these methods are used for stiff problems, 
they suffer from an order reduction phenomenon. Specifically, if p is the largest integer for which 
the following holds: 


( 1 . 1 ) 

and: 


( 1 . 2 ) 


t> 


y ‘(0 

m 


T[e; Ae ; . . . ; A p ~ 2 e\ = [Ae; 2A 2 e; . . . ; (p - l)A p-1 e] 

then it often happens that only a fc nun (* /, p) type convergence can be proved or demonstrated com- 
putationally. (See e. g., [6], [7], [4], [10], [12], [15], and [11]) Also, note the barrier: p < q + 1. 
With regard to stability, let r(z) be a rational approximation to the exponential e~ z defined by: 

(1.3) r(z) = 1 — zb T (I + zA)~ 1 e. 
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An IRKM is said to be Ao-stable if: 



kWI < i 

Vz > 0, 

/-stable if: 

I'M | < 1 

V5R{*} = 0, 

A-stable if: 

K*)I<1 

VSft{z} > 0, 

and algebraically stable if: 




B and M are positive semidefinite. 


Note that the last notion of stability is the strongest among these. [7] In fact, for the methods 
which are algebraically stable with B actually positive, there exist results for certain parabolic 
equations which guarantee decay of approximations with respect to the time step. [10] 

Now, concerning implementability, there has been much effort devoted to the development of 
IRKM’s for which the eigenvalues of A are identical and real. ([8], [9]) As indicated in the next 
section, these so-called singly implicit (SIRK) methods offer a computational advantage over other 
IRKM’s on serial machines. However, it is explained in section 2 that in a parallel environment the 
preferred methods are those for which the eigenvalues of A are distinct. In this work, the latter are 
referred to as multiply implicit (MIRK) methods. Further, they are called real if cr(A) C R, and 
otherwise complex. 

It can be seen in section 2, that real MIRK’s permit a greater degree of parallelism than those 
which are complex. However, if the eigenvalues of A are real, then v < q + 1. [14] Nevertheless, 
while the classical order can be increased at the cost of introducing complex eigenvalues, the order 
reduction phenomenon enforces the p < q -+■ 1 barrier for the problems which motivate this work. 
Hence, the principal interest here is in real MIRK’s. 

In this connection, it is important for certain parabolic problems that for every positive integer 
q, there exists a g-stage, real MIRK which is Ao-stable with maximal order v = q+ 1 and p — q + 1. 
This fact follows from results in [2],* but an independent and direct proof is given in section 4. 
Also, for hyperbolic problems, it is shown in section 5 that for every even integer q, there is a 
g-stage, real MIRK which is /-stable with order v — q. For a related result, see Bales, Karakashian, 
and Serbin.* 

Concerning methods which are more stable, Wanner, Hairer, and Nprsett [16] have established 
that if an A-stable SIRK has order v — q + 1, then q = 1,2, 3, or 5. Also, for SIRK’s which are 
actually algebraically stable, the limit p = q + 1 is achievable only for q < 2. [3] Nevertheless, it is 
shown in section 3, that for a wide range of problems, an IRKM can be modified in such a way that 
the order reduction is no worse than even if p < g + 1. Hence, without regard for the 

value of p, one is led to ask about the existence of real MIRK’s which are A-stable with maximal 
order u = q+ 1. In section 4, it is shown that a necessary condition on such methods is that q = 1, 
2, 3, or 5, i.e., the result of [16] is generalized to the case of real, distinct eigenvalues. Then in 
section 6, some useful examples of algebraically stable MIRK’s are presented. 


‘Bales, Karakashian, and Serbin, private communication. 

* BALES, L. A., KARAKASHIAN, O. A., SERBIN, S. M., On the Stability of Rational Approximations to the 
Cosine with Only Real Poles. (To appear.) 
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2 Parallel Implementation. 

The primary purpose of this section is to support the claim that among IRKM’s, MIRK’s 
are the preferred methods in a parallel computing environment. For definiteness, let S h be a 
finite dimensional function space and suppose that an approximation is required for the solution 
u : [0, i*] — ► Sh, to the initial value problem: 

f Dtu = — Lh(t)u + f(t, u) 0 < t < t* 

\ u(0) = u° 

where h amounts to a stiffness parameter, and for simplicity, f(t, u) is assumed to be smooth and 
independent of h. Also Lh{t) is assumed to be linear and selfadjoint with positive spectrum. Such a 
problem could of course arise from the semidiscretization of a semilinear parabolic initial boundary 
value problem. For its temporal discretization, let t* = kn * , t n = kn, and make the following 
definitions: 

LI = L h (t n ), Zl = diag{L£}, ft = diag {L h (t n + kn)}. 

qX9 l<t<g 

It is well-known that implementing (1.2) as it stands can be very expensive because of the burden 
involved in computing the stages. Nevertheless, for simplicity here, assume that for the first v — 1 
time steps, the stages are indeed computed exactly. Then, moving toward a cheaper procedure, 
given {/(t n_m , U n ~ m )y^} Q , define approximate stages U n = (U n ' 1 ,U n ' 2 , . . . , U n,q ) T according to: 

U n = [/+ kAZl\~ l {eU n + kAE n f} 


where: 


m=0 


v — 1 

^ a im m' = (-l) i ejT J e 0 < / < v - 1, 1 < j < q, {m l \ m=t = Q = l), 

m= 0 

and terms such as A£% for example, are understood in the sense of composition of operators defined 
on [Sh] q - Even though this extrapolation circumvents the need for solving a system of nonlinear 
algebraic equations at every time step, inverting the full operator [I + kAZ^} could still be very 
costly. Hence, suppose: 

A=S~ 1 \S, 

A = diag {A,} + subdiag {0,} A,- > 0, 1 < t < q, 6i = 0 or 1, 2 < i < q 

1 <i<9 2<»<g 

so that Vj” « U n can be obtained by the (outer) iterations: 


[/ + fcA£^](5Vi n ) = {SeU n + kSA(Z% - Z^V^ + kSAE n f} = R? 1 < l < l n 


where: 


min K‘0 / . \ 

Vq” = / nun\n,V) j 

m= 1 V m ) 


1 < n < n* — 1, 


V 0 ° = eU° 


= eu 


0 


> 
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and arg com p U ^ ec j as indicated below. Now, consider the simple but important 

observation that if A,- ^ A j, i ^ j and $i = 0, 2 < t < q, then for arbitrarily large q, the block 
system above completely decouples into the following equations which can be solved in parallel: 

[/ + fcAiLEKSV?*),- = (*,"), ■ !<«'<?• 


On the other hand, if a SIRK is used for which A,- = A, 1 < t < q and 0 t - = 1, 2 < t < q, then at each 
time step, the approximate solution is required of systems involving only a single new coefficient 
matrix with the dimension d h of S^. However, Karakashian and Rust * have compared a two-stage 
MIRK and a two-stage SIRK on a two-processor IBM-3081D at the University of Tennessee, and a 
four-processor CRAY-XMP48 at the super computing center in Pittsburgh. This work suggests that 
as the dimension dh increases, making the solution of the time stepping equations the dominant 
operation, the speed-up quickly exceeds unity and then approaches the ideal factor of q. Therefore, 
if a certain sufficiently large dh (high spatial accuracy) is not required for a given problem, it 
could of course happen that the serial calculation would offer superior performance for a fixed q. 
Nevertheless, the temporal accuracy requirement could be satisfied by taking q large enough that 
a dramatic reduction in execution time is achieved with parallel computations. 

Finally, in any of these cases, the factorization of new coefficient matrices at every time step can 
be avoided by using a preconditioned iterative method to approximate Vj n with (inner) iterates, 
say {kj”}o<j<;„- Further, for certain problems, it can be shown [10] that there exist integers m n 
and j n such that: 

1 n* — X 

~7 ]C m njn <C ^ C(h, k ) 
n n=0 

while the convergence order obtained for the scheme: 


f [P+ 1 
{ U n 


(1 - b T A~ 1 e)U n + kb T A~ 1 U n 


tJn 


0 < n < n* — 1 


is the same as that obtained by using any additional outer or inner iterations. 

Studying methods similar to those discussed above, certain authors have proved convergence 
estimates for partial and stiff ordinary differential equations obtaining 0(fc mm ( l/,p )) for the temporal 
discretization. (See e. g., [6], [10], and [4].) However, as advertised in the Introduction, for a 
wide range of evolution equations, there is a systematic way of modifying an IRKM so that order 
reduction can be avoided. Specifically, it is shown in [12] that for semilinear parabolic equations 
with time independent coefficients, order reduction can be completely eliminated by computing the 
stages using an extrapolation procedure in which the constants are defined by: 


v— 1 

= (-1 ) l l\ejA l e 0 < / < v - 1, 1 < j < q. 

m=0 

It is also shown in [12] that for linear parabolic equations with time dependent coefficients, the 
order reduction can be rendered no worse that 0 (A: min ( l ' ,?+1 )) by computing the stages as indicated 

‘KARAKASHIAN, 0. A., Rust, W., On the Parallel Implementation of Implicit Runge-Kutta Methods. (To 
Appear.) 
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at the beginning of this section but with ££ replaced by: 


min(i/— 1,9) 

E 

m=0 


2; = E r m £;- m 


where: 


and: 


min(i/— l,g) 

£ r m m' = (-1 ) l D l 

m=0 


0 < / < min(i/ — 1, q) 


D[e\ Ae; ...;A q 1 e] = [Ae\ 2 A 2 e ; . . . ; qA 9 e]. 

Finally, in [ 12 ] it is shown that an approach similar to this can be used for quasilinear problems to 
obtain comparable results. Note that in each case, a special starting scheme is used to initiate the 
indicated extrapolation but the details are not provided here. 

Now consider the case that <r(A) (£. R, but assume that A can be transformed to quasidiagonal 
form: 

A=S' _1 AS, A = diag{A<}, 1 < m < <7 

l<i<m 


and either: 


A,- — A,- > 0 or A,- = 


A 


a, 

A <*» 


> A - > 0 


1 < 1 < m. 


Then, the block linear systems discussed above decouple to equations of the form: 

[l + k\LM = 4> 

and: 


’ I + kocLl 

1 

V’l 


<f > 1 

mi 

I + kaL% J 

. ^ 2 . 


<t> 2 


Again the subordinate equations can in principle, be solved simultaneously. Further, the solutions 
ip 1 and tp2 for the 2x2 system can be computed in parallel according to: 

[I + ZakLl + (a 2 + p 2 ){kLl) 2 }rPx = [/ + kaLfifr. + kfiLlfa = X i 

and: 

[I + 2 akL% + (a 2 + ^){kL n h ) 2 }rP 2 = [I + kocL n h \<h - kf}Ll<h = X2- 
Following Baker, Bramble and Thomee [l], since for x € R: 


= SR 


zi 


Ll + z 2 x 


z\ = 1 - a /3 1 », z 2 = a + ( 3 i, 


1 + 2 ax + (a 2 + /? 2 )x 2 
complex arithmetic can be used to obtain: 

t pi = + zikLhY'xi) 


* = 1 , 2 . 


Note that as mentioned in the Introduction, the depth of decoupling which the real MIRK’s 
allow is greater than their complex counterparts. Nevertheless, examples of each are given in 
sections 4 r 6 . 
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3 A Barrier for A-S table Real MIRK’s with Maximal Order. 


In this section, it is proved that a necessary condition for a g-stage, real MIRK to be A-stable 
with maximal order v = q + 1, is that q = I, 2, 3, or 5. The plan of the proof is roughly as 
follows. According to the work of Nprsett and Wanner [13], for a real MIRK to be A-stable with 
maximal order, the eigenvalues of A must be linked to a certain type of hypersurface near the 
middle of the so-called real-pole sandwich. However, it is shown below that when constrained to 
such hypersurfaces, the rational function (1.3) has absolute value bounded by unity for infinite 
argument only if q = 1, 2, 3, or 5. The latter claim is established by considering this absolute 
value as a function depending on (the reciprocals of) the eigenvalues of A, and showing that it is 
minimized over (the positive part of) a maximal order hypersurface of the real-pole sandwich at 
the point where the eigenvalues are the same. Then since such points have been studied carefully 
for SIRK’s, the remainder of the proof follows from the work of Wanner, Hairer, and N0rsett. [16] 
Before proceeding with the details, the appropriate notation is developed. First, define the 
so-called symmetric polynomials for x G R m : 

5 0 (x) = 1, S/(x) = £ x h • x «i 1 < f < m > ■S’ m+1 (x) = 0. 

Next, some relevant properties of these polynomials are enumerated. By direct calculation: 

S/(x) = Xi £ x h x u - 1 + £ x h ®«i 

(3.1) 


Thus: 


Xid x .Si(pc ) + d Xi S t+ i(x) 


1 < t < m, 0 < / < m Vx e R m . 


(3.2) d x .Si + i{x) = ^(xJIjj.^o 1 < » < m, 0</<m 

Also with e m € R m having all unit coordinates: 


Vx e R m . 


Si(xe m ) = x l £ l = x ‘(7) 

*1 ' / 


(3.3) 

since the sum is the number of combinations of m distinct integers taken / at a time. Finally with 
x' 1 = (xi 1 ,X 2 1 ,..., x ~ 1 ) I ' : 

Sl( x ) 


(3.4) 


| S'm( x ) 


= S m -l{x X ) 


0 < l< m 


Vx -1 e n m . 


Now according to [13] and (3.4), the rational function (1.3) has real poles and order of approx- 
imation to e~ z equal to at least q,* if and only if it has the form: 

9 l 1 ilm 1 


(3.5) 


r(z) = 


n-‘Y t ii-rt ^ 

i=0 m= 0 V 7, V- i =0 f7i=0'' >' 


T-Sii*- 1 )*' 

l=o 


£s,_,(x)2' 

1=0 


"This notion of order is of course weaker than that described in the Introduction for IRKM’s, but the distinction 
should be clear from the context. 
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for some x 1 £ R 9 . Also since: 


E s h( x ) 2 '= II( x »" +2: ) 

1=0 m= 1 

the poles of r(z) are So it follows from work in [13] that a g-stage, real MIRK can be 

A-stable with maximal order q + 1 only if x is contained in one of the q hypersurfaces of the set;t 

W = {xe R*: Xi > 0, = 

For convenience, these sheets or connected components of W are ordered as follows. First, define 
the Laguerre polynomial of degree m: 

= E(-i>'(7)£ 


and the values by: 

D = o, 


0 < x™ < x™ < 


< 


x 


m 

m* 


Then using (3.3): 


W = Wi U W 2 U . . . U W q with exf € Wi, 1 < i < q. 

Now, the following is an adaptation of Theorem 10 of [13]. 

Lemma 3.1 Let an A-atable, q-stage IRKM be given with a rational function ( 1.8 ) having only 
real poles and maximal order q + 1. Then (IS) must have the form (8.5) and it is necessary that: 


X £ W q + 1 U W q - 1 if q is odd , or 

2 2 


x £ Wg if q is even 
2 


where Wo = 0. ■ 

The next lemma is proved in section 5 of [16]. 


Lemma 3.2 For m, i > 1, | L m (x™)\ < 1 only if: 


m > { 


1, 

5, 

9, 

6 * - 10 , 


* = 1 
i = 2 

i'=3 
* > 4. 


E 


*The maximal order hypersurfaces of the real-pole sandwich are defined in [13] as the sheets of M = {y £ R ? : 
l-o = 0); so by (3.4), to every point in W there corresponds a point in M with reciprocal coordinates. 
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Recalling (3.5), define: 


*(x) = = r (°°)- 

1=0 

Next, let an IRKM be given as specified in Lemma 3.1, so that in particular, ]r(oo)| < 1. Now 
suppose that the following has been established: 

(3.6) ^rnf |£(x)| = \L q {x1)\ 1 <i<q. 

Then according to Lemma 3.1, one of the following must hold: 


q is odd, 

xe W q + 1 , 
2 

and: 

\Lg(x q q+1 )| 
2 

< 


q is odd, 

XG Wg-l, 
2 

and: 

\L q {^)\ 

2 

< 

> l^(x)l < 1 

q is even 

xe W ± , 

and: 

\L q (x\)\ 

< 



2 2 


Finally according to Lemma 3.2, the first case implies that q — 1, the second that q = 3 or 5, and 
the third that q — 2. This is the advertised result and what remains is to establish (3.6). 

Lemma 3.3 The constrained critical points of R(x) on W are precisely those points in W with 
identical coordinates. 


Proof. First, form the Lagrangian: 

F(x, A) = £ m^$(x) - A g 
so that if £(x) has a constrained critical point at x* 6 W, then: 

(3.7) o = a Ii F(x%A) = -£i^a,,.s, +1 (x-) + A£iJ| T a 1 .s, +1 (x-) i <•'<?. 


Also by (3.1): 


o = E^(x-) = 


(3.8) 


= 


9 - (-1)' 


9-1 (-i y 
(/+ 1) 

9-1 (-iy 


1 < i < q. 


Now for the sake of contradiction suppose that for some j: 

(»•») = °- 


8 



Then it follows from (3.8) that in addition: 


(3.10) 


j),d Iy S i+1 (x*) -0. 


Using (3.2) and (3.4), it follows that if the coordinates of y* G R 9 1 are reciprocal to those of x* 
other that a;*., then: 

My- 1 ) „ ... 

3.,S,(x>) S q . 1 (y- 1 ) ’■ 

So, dividing (3.9) and (3.10) by d Xj .S q (x*) gives the following conditions: 


9-1 r_iV 

E 77-T^T S «-l-<(y') = 0 


J 0 ((+l)! V, - ,(y) - 


y G R" 


However, as indicated in [13], with: 

(3.11) No(t\y) = ^(— l) n-z 5 n _j(y)^, N m+ i(t;y) = N m {s;y)d. 

the following must be satisfied: 

{yGR n :iV 1 (l;y) = 0} n {y G R n : JV 2 (l;y) = 0} = 0. 

Hence, (3.9) and (3.10) cannot hold simultaneously. Therefore, by (3.7) and (3.8), = A = x*-, 

1 <i,j<q- ■ 

Now in order to capture the global minima of |£(x)| on the components of W, the following 
corollary is concerned with constrained critical points on the boundary of W . Its proof involves 
only appropriate adjustments in the dimensions of the above argument. 

Corollary 3.1 Given a multi-index m = (mi, m 2 , . . . , m q ) where m,- = 0 or 1, 1 < 1 < q, and 
|m| = ™i> 1 < |m| < q, define: 

R™ = {x G R 9 : Xi = 0, if mi = 0, Xj > 0, if mj = 1, 1 < i,j < q}. 

Then the constrained critical points of £(x) on W PlR^ 1 are precisely those points in W PiR™ with 
identical coordinates. m 

The next lemma is not actually contained in Theorem 8 of [13], but the argument required is 
essentially the same. Nevertheless, a brief proof is offered for completeness. 

Lemma 3.4 Let m = |m|, 1 < m < q, and suppose e m G R™ has only unit nonzero coordinates. 
Then the points ofWn R™ with identical coordinates are precisely: 


x m e m € W; 


1 < 1 < m. 


and: 


Z(x?e m ) = L m (x?) 


1 < » < m. 
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Proof. Fix m, m = |m|, with 2 < m < q, and n, n = |n|, with n = m — 1. Also, suppose 
R“ C R™ and that e m G R™ and e 11 G R+ have only unit nonzero coordinates. Next, with 
A(0) = {x G R™ : x = xe m ,x > 0} and -X’(l) = {x G R+ : x = xe n ,x > 0}, suppose that 
for 6 G [0, 1), X(0) is a ray from the origin into R™ passing smoothly from -X’(O) to A(l). With 
(3.3), it follows that *(0) D W = and X(l) n W = {a:r ell }"=i- So for 6 e [0,1), let 

{x,(0)}J^ 1 denote the points of X(0) n W ordered according to increasing magnitude. That these 
points remain separated, preserved in number and order, and bounded for 9 in compact subintervals 
of [0,1), can be seen by noting that the explicit calculation of their coordinates involves the roots 
of a polynomial which has degree m while 9 < 1. Since the components of W never intersect, [13] 
these roots can never coalesce to become multiple or complex or to change relative positions. Now, 
since their number decreases by one as 0 passes to 1, it only remains to determine whether x m 
escapes to oo or Xi passes through the origin. Since W n 0 = 0, the indicated alignment between 
the points and {W,-}i<,-< 9 is established. Finally, the values for £(x) are obtained 

with (3.3). ■ 

Now the following shows that the constrained critical values of £(x) determine the constrained 
global minima of |£(x)|. 

Lemma 3.5 The constrained global minima of | Z (x) | on the components ofW are given as follows: 

inf |£( x )l = min {|L m (*P)|} 1 < * < ?• 

xG^,- %<m<q 

Proof First, note that £(x) never vanishes on W , because the rational function (3.5) cannot have 
degree q — 1 in the numerator and order q+1 simultaneously. [14] Hence, on each W i nR+S |S(x)| 
is a smooth, positive polynomial with exactly one critical point where it must be minimized. The 
rest of the result follows with Corollary 3.1 and Lemma 3.4. ■ 

Now, the following lemma yields (3.6). 

Lemma 3.6 The Laguerre polynomials {L m (x)} m >i, and the values satisfy: 

(3.12) IMxDI < lim-x^r 1 )! l<*<m — 1, m > 2. 

Proof. By the identity: 

(3.13) Lm(x) = L m +i(x) — m --j- L m+ 1 (x) 
the inequalities of (3.12) are equivalent to: 

|L m+1 (xDI < |im(xr _1 )| 1 < i < m - 1, m > 2 

which can be established by showing that: 

(3.14) L m (x7 l_1 ) < L m (xf l ) = L m+ i(xD < 0 for odd t < m - 1 
and: 

(3.15) 0 < L m+ i(xf l ) = L m (x- n ) < L m (xJ* _1 ) for even i < m - 1. 
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First note that since L m+ i(x) has exactly m + 1 simple real roots, it has exactly m local extrema 
which must be associated with the m simple real roots of L' m+1 (x). Therefore since L m +i(0) = 1: 


(3.16) (-l)*'L m+1 (*n > 0, 

So using this with (3.13): 


(-i)‘ +, c + i(*n > o 


1 < » < m. 


V +1 T m r m 

(-lyiUxT) = >0 1 < t < m. 


i«+l x m 


Now assume that i < m — 1 is odd. Then on the interval [xj 71 , X™ j], £»»(*) is decreasing at the 
beginning and increasing at the end, while L' m (x) vanishes only at x™~ x € (if 1 ) 1 ™!) where L m (x) 
is minimized. Combining this fact with (3.13) and (3.16) gives (3.14). Also (3.15) follows similarly. 


Finally, the result of this section is summarized in the following. 

Theorem 3.1 Let an A-stable, q-stage IRKM be given with a rational function (1.8) having only 
real poles and maximal order q + 1. Then q = 1, 2, 3, or 5. ■ 


4 A Family of High Order ylo-Stable Real MIRK’s. 


In spite of the barrier established above, a family of IRKM’s is constructed in this section, 
which contains for every positive integer q, a g-stage, real MIRK which is Ao-stable with maximal 
order q + 1. Such methods are useful for parabolic equations, since the following is pivotal in the 
stability analysis: 

MM*)II < 1 

where || • || is an appropriate norm. [10] Note that if Lh[t) is selfadjoint and positive definite, 
then this inequality follows from Ao-stability and a spectral argument. Otherwise, a restrictive 
relationship between h and k must be imposed. 

Toward the goal of this section, recall the functions {iV m (t; y)} m >o defined in (3.11). Then let 
e= (ei, e-i , . . . , e q ) T be a unit vector chosen so that: 


e,- > 0 and e,- ^ ey 1 < t ^ j < q. 

According to [13], there are exactly q positive solutions to: 

Ni(l,ye) = 0. 


So if y* is the largest, set y* = y*e. Then it follows that the function: 

Ni(t,f) _ t^A^M-V) 

A((0,y*) (-l)«S,(y*) 

vanishes at t = 0 and t = 1 but for no t G (0,1). Further, since G'( 0) = 1, it follows that for every 
t € [0,1], t?"(t) < 0. Finally, note that by (3.4): 



l=o 


(-i) 1 s q -i(r) 

" S q ( y*) 


sty*- 1 )- 
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Lemma 4.1 |G'(1)| < 1. 

Proof. As explained in the proof of Lemma 3.5, £(x) never vanishes on W . Also, from the proof 
of Lemma 3.4, it can be seen that W\ is compact. So, the global maximum of |£(x)| on Wi is 
determined by the constrained critical values of £(x) given by Corollary 3.1 and Lemma 3.4. Hence 
by Lemma 3.6: 

sup |£(x)| = max {[Lm^DI} = \ L i{ x i)\- 


xetVi 


l<m<g 


Then since some calculations show that y* 1 E Wi, |G !/ (1) | < |Li(xJ)| = y/2 — 1. 

Now set: 

ti-ft 

R{z) = i=2 ^ - 

y*) zl 

1=0 

taking Q(z) as the denominator. 

Theorem 4.1 Any IRKM for which (1.8) is equal to R(z) above, must be Ao-stable. 
Proof. First, note the error formula of [13]: 


R(z) = 

With x > 0, integration by parts gives 





f-xl ? 

= i+ Q ( x) r 


. , r 

Q(x) 

or: 

R(x) = 

where: 

0 = f Sqir) 6 [0,1), 


^Z s i{ y*) x 

1=0 


E [0, 1), and J = [\l - e~ x ^)G"{t)dt. 
Jo 


As indicated above, G''(t) < 0 and |G'(1)| < 1, so: 

-J = \ J\ < - f 1 G' f (t)dt = 1 - G'(l) < 2. 
Jo 


Hence: 


|i?(x)| = R{x) = e~ x (l - 0) + 0(1 + J) < 1 


Vx > 0. 


Now the next three theorems show that R(z) can be used to construct a g-stage, real MIRK 
which is Ao-stable with order q + 1. The first is due to Bales, Karakashian, and Serbin [2]. 
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Theorem 4.2 If the coordinates of y* £ R ? are distinct and nonzero, then No[t\y*) has q distinct 
real roots . ■ 

Theorem 4.3 Given distinct real roots {r,}? =1 of No(t;y*), the following are well-defined: 


A = TVRV - 1 and b = {V^Re. 

and the resulting IRKM has order at least q. Furthermore, if N\(l,y*) — 0, then the order is q + 1. 
Proof. With the conditions of Butcher: 

B(0 : lb T T l ~ l e = 1 1 < l < £, 

C(f) : lAT l ~ l e = T*e 1 < / < f, 

conditions B(u) and C(/i) imply order v if v < fi+ 1. [5] Now the first statement follows since the 
assumptions amount to conditions C(q) and B(q). Then if ^(^y*) = 0, by B(q): 

= «>».(», y*) - 

1=0 i= 1 

= E 6 ‘W-?^o(r,,y*)] = b T T<*e 

»= l 

and the second statement follows with B(q + 1) and C(q). m 

Finally, the next theorem follows from Theorems 2 and 4 of [13]. 

Theorem 4.4 Given an IRKM constructed as described in Theorem 4-8, the rational function (1-8) 
must have the form (8.5) with x -1 =y* . ■ 

5 A Family of High Order I - Stable Real MIRK’s. 

In this section, a family of IRKM’s is constructed which contains for every even q, a g-stage, 
real MIRK which is /-stable with order q. As briefly indicated below, such methods are useful for 
hyperbolic problems. 

Suppose that an approximation is required for the solution ti : [0, t*] — ♦ Sh, to the initial value 
problem: 

D\u = -L h (t)u+ f(t,u) 0 <t<t* 

u(0) = u° 

D t u( 0) = u 1 
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where Lh{t) and f(t,u) are assumed to have the same properties as mentioned in section 2. Note 
that this problem can be expressed in first order form as follows: 



As opposed to the parabolic case, with: 



the pivotal stability inequality here is: 

IIK*L»||| < 1 


where ||| • ||| is an appropriate norm. Since L" has only imaginary eigenvalues, a certain spectral 
argument shows that /-stability is crucial. [11] 

Now in the sequel, q is implicitly assumed to be even. Toward the goal of this section, recall the 
functions y)} m >o defined in (3.11). Then let e = (ei,ei,. . . ,e q ) T be a unit vector chosen so 

that: 

e,- = — e . £ > 0 and e,- ^ e ; 1 < t j < |. 

According to [13], there are q real solutions to: 

Ao(l,ye) = 0 

and exactly £ of them are positive. So if yo is the largest, choose y* > yo and set y* = y*e. Then 
it follows that: 

rm = Mil!) - ^o(M'V) 

U N o (0,y*) (-l)§|S,(y*)| 

decreases in a strictly monotonic fashion from 1 to G( 1) > 0 as t passes from 0 to 1. Now set: 


<3M= n [1 - (sM ! l = I>(y V, 

m=l 1=0 


9 1 

p(.) = £(-*)' J2 ±-L-s m ( y -) 

1=0 m= 0 V >' 


and take R(z) = P(z)/Q(z). 

Theorem 5.1 Any IRKM for which (IS) is equal to R(z) above, must be I-stable. 
Proof. First, recall the error formula of [13]: 


R{z) = e 




So with z = iy: 


Wv)l = 


l-ty 


n 


( yy * m ) 2 


i [i + (yy* m y 


f 1 e - ivi G(t)dt 

Jo 
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Integration by parts gives: 


|i?(iy)|=|(l-0) + 5(e-‘'G(l)-J)| 


where: 




By construction, G'(t) < 0 , so: 


|J| < - f l G\t)dt = 1 -G(l). 

Jo 


Further, since G( 1) > 0: 


|fl(z)| <(!-*) + 9G{ 1) + 0(1 - G(l)) < 1. 


Now, in view of Theorems 4.2 - 4.4, for every even integer q, the rational function R(z ) above 
can be used to construct a 9 -stage, real MIRK which is /-stable with order q. 

6 Examples of Algebraically Stable MIRK’s. 

In this section, certain methods are presented which have the properties discussed in previous 
sections. In the discussion, use is made of Theorem 6.1 below, which is due to Hairer and Wanner 
[9]. For this, some relevant constructions are now provided. Given B > 0, let C be determined by 
the Cholesky decomposition: 


Then take: 


so that: 


Finally with: 


V t BV - C T C. 


W = VC' 


W t BW = I. 


{ *’ ' 
( 0, o 


1 < m 
otherwise 


let I m = H = {(* + J “ 1 ) 1 }?.,=!> and: 


X G = diag{|, 0, 0, . . . , 0} + subdiag{^,} - supdiag{^,}, 

qxq 1 <»< ij -1 

Now the following can be stated. 

Theorem 6.1 Let an IRKM be given for which B > 0 and: 

IiV T BVIj = IiHIj i + j-1 
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Let m = | — 1)] and suppose X is a q X q matrix satisfying: 

XI m = X G I m , I m X = I m X G 


and if v is even: 

^Tn+iXem+i = 0 . 


Then if A= WXW~ l the IRKM has order v. 


First, the following family of two-stage methods has been described by Karakashian [10]: 


ti{t 2 - |ri) 

_l r 2 

2 r l 


r 2 - ti 

T2 ~ Tl 

•1 

ir 2 

2 r 2 

T-2(^2 - Tl) 

Ti 

7- 2 - Ti 

T 2 - Tl 

t 2 -\ 



T 2 - Tl 

Ti ~ Tl 



The W -transformation leads to: 


n = (Ai + A 2 ) — (Aj + A 2 ) 2 
t 2 = (Ai + A 2 ) + (Af + Aj ) 2 


Ai > 2 > 


A2 = 


Al— ; 



' 1 2v/3(t-i - 

*)' 


x 1-1 

2 - 2Vf 

w = 



and X = 



1 2v/3(r 2 - 

1). 


1 

L 273 1 


where: 

x = Ai + A2 — j > 0. 


Since B > 0 and V T BV = H, by Theorem 6.1, these methods have maximal order v = 3. In 
fact, p = 3 as well. With regard to stability, since W T BW — I, the algebraic stability matrix M 
satisfies: 


W t MW = W T BWW~ 1 AW + (W~ 1 AW) T W T BW + W T BW[W- 1 e(W- 1 e) T ]W T BW 


= X + X T — eiej = diag{0,2i}. 

Hence, these methods are algebraically stable. Also in connection with section 2, note that A has 
real, distinct eigenvalues: 


A = S~ 1 hS, 



s-' = 


Tl — A 2 Tl — Ai 

J"2 — A2 7"2 Ai 


Next, the following family of three-stage methods has been described in [12]. First, choose Ai, 
A2, and A3 distinctly and to satisfy: 


Ai, A 2 , A 3 > 0 Ai + A 2 + A 3 — J- > 0 and A 3 


j-AiA 2 - |(Ai + A 2 ) + ^ 
AiA 2 — |(Ai + A2) + | 
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For example, it is sufficient to choose A3 as indicated after taking: 

1 \ _ 1 

Ai > | and A2 > - — - — p- 

Ai — j 

Now define: 


xi = Ai + A2 + A3 — |, x 2 = [(A1A2 + Ai A3 + A2A3) — j(Ai + A2 + A3) + |] 2 , 


and: 

a,- = A,- - xi, = j - A,- 1 < t < 3 . 

Also take A = diag {A,}, 



1 

2 

_ 1 
2%/3 

0 


01 1 

a 2 

a 3 

x = 

1 

2\/3 

0 

-x 2 

, and Y — 

2 \/ 3 q!i/?i 

2 V*cc 2 ( 3 2 

2 \Z%Ol 3 p 3 


0 

*2 

Xl 


2s/2x 2 Pi 

2 \fZx 2 {i 2 

I 1 

to 

CO 

1 , , 


Then X = Y AY Now choose n, r2, and r 3 to satisfy: 

( r l ~ r i + k)( r 2 ~ r 2 + g) , . _ frr 2 - |(r x + r 2 ) + \ 

7- — < 0 and T3 — - f-. 

T l T 2 ~ 2 ( r l + T i) + 3 t 1 t 2 - |(ri + r 2 ) + i 

For example, it is sufficient to choose T3 as indicated after taking: 


1 a j 1/, 1 a 

ri< 2 ( 1 “vS> “ d 2 (1+ ^ )<rj ' 


Then it can be shown that: 

7 

a = \t x t 2 T3 - \(t x t 2 + t x t 3 + t 2 t 3 ) + i ( r x + t 2 + r 3 ) > — 


36 


so take r = (a — 2 and define W by: 


1 2\/Z{t x -\) (rf - ri + |)/r 

W = I 1 2 \/ 3 (f* - £) (r| - r 2 + |)/r 

[ 1 2 v/ 3 (r 3 - |) (rf - 7-3 + |)/r 

Next, since it can be shown that ry 7^ ry, 1 7^ j, choose the vector 6 according to: 

b = (V T ) -1 i?e. 

Finally take: 


A = \VXW ~ 1 = S~ 1 A.S, 


S = Y-'W - 1 
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and note that A has real, distinct eigenvalues. Concerning stability, note that B > 0 since V T BV = 
H + diag{0, 0, a — g-}. Also since W T BW = /, the algebraic stability matrix M satisfies: 

W t MW = W t BWW~ 1 AW + (W~ 1 AW) T W T BW + W T BW[W~ 1 e{W- 1 e) T }W T BW 
= X + X T — eief = diag{0, 0, 2ii}. 

Hence these methods are algebraically stable. Now since B > 0 and V T BV is as indicated above, 
it follows from Theorem 6.1 that these methods have maximal order v = 4. However, p = 2. 
Nevertheless, it is now shown that for each method in the family, there exists a matrix D as defined 
in section 2: 

£>[e;Ae;A 2 e] = [Ae; 2A 2 e; 3A s e] 

which actually eliminates order reduction since u = 4 = q + 1. For this, suppose that there are 
constants such that: 

cie + C2Ae + C3A 2 e = 0. 

Since a fourth order method satisfies: 

l\b T A l_1 e = 1 1 < / < 4 

it follows that: 

ci + |c 2 + |c 3 = 0 

\ci + \c 2 + ^c 3 = 0. 

Therefore 12ci = — 2c 2 = c 3 , and if these constants are nontrivial: 

[/ - 6A + 12A 2 ]Se = 0. 

However, since 1 - 6t + 12t 2 > 0, Vt £ R, it is necessarily the case that ci = c 2 = c 3 = 0. Thus the 
matrix [e; Ae; A 2 e] is invertible. 

It is also possible to construct similar methods which have five stages and order six, but the 
details are not provided here. Instead, some useful calculations for some well-known complex 
MIRK’s are given. Recall from section 2 that at least for semilinear problems, order reduction 
can be eliminated; so for such problems, there is good reason to consider methods of the following 
form. First, the Gauss-Legendre methods are algebraically stable with v = 2 q. From this family, 
the three-stage method follows. 


5 

80 - 24\/l5 

50 - 12>/T5 

5 - \/l5 

36 

360 

360 

10 

50+ 15\/l5 

2 

50 - 15>/l5 

1 

360 

9 

360 

2 

50+ 12\/l5 

80 + 24\/l5 

5 

5 + y/lE 

360 

360 

36 

10 

5 

4 

5 


18 

9 

18 
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Also, the ^-transformation leads to: 



' 1 

3 

2 1 
75 


r l 
2 

_ 1 
2V3 

0 

w = 

1 

0 

2 

and X = 

1 

2\/3 

0 

_ 1 
2s/lt 


1 

3 

2 

7E J 


0 

1 

2v^5 

0 


In connection with section 2, note that: 

A = S~ 1 AS S^Y-'W- 1 


where X = FAF -1 , 



' A 

0 0 


2>/l5A 

-12a/? 

12a 2 — 6a + \ 

A = 

0 

a —/3 

, Y = 

12v/5A(i- A) 

~\/3/? 

— y/Za 


0 

/? a 


2v/3(|-A) 

0 

l 

2-s/5 


and the above constants are given by: 

A = ^-[(5 + 3^5)? + (5 - 3v/5)?] + i 



20 ? 

120 


[(5 + 3v/5)? + (5-3v/5)?], 


and /? = v/3y^[(5 + 3-\/5)? - (5 - 3\/5)3] 


Now the Gauss-Legendre methods fail to satisfy |r (oo)j < 1, and the latter condition is required 
in the analysis of [12] for certain nonlinear problems. On the other hand for example, the Radau 
IIA methods satisfy this constraint. Also, they are algebraically stable with v = 2q — 1. From this 
family, the three-stage method follows. 


88- 7\/6 

296 - 169\/6 

-2 + 3\/6 

4->/6 

360 

1800 

225 

10 

296 + 169\/6 

88 + 7\/6 

-2 - 3\/6 

4 + Ve 

1800 

360 

225 

10 

16- yf& 

16 + \/6 

1 

l 

36 

36 

9 


16- y/6 

16 + \/6 

1 


36 

36 

9 
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Also, the W -transformation leads to: 



1 

*51 
00 
1 in 

■s 

i 

— 2\/5+3'\/30 
25 


r 1 
2 

_ 1 
2V3 

0 

w = 

1 

-■\/3+3^ 

5 

— 2\/5— 3\/30 
25 

and X = 

1 

2\/3 

0 

1 

2\7l5 


1 


V/5 


0 

1 

2\/l5 

1 

10 


In connection with section 2, note that: 

a = s -1 as s = y _1 py _1 

where X = YAY -1 , 



' A 

0 0 


2VI5(A - i) 

0 

-V5 

A = 

0 

a — / 8 

, Y = 

12\/5(A — j^)(| — A) 

2/?\/l5 

(2a - l)\/l5 


. 0 

(3 a 


2^3(1 -A) 

1 

•Hi'* 

o 

l + 120(a-i)(a-i) . 


and the above constants are given by: 

A s i(3t - 3-») + I Q = I-i( 3 r-3-i), and ^=l(3f + 3 i). 
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